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Abstract 

We study the thermal conductivity of the one dimensional Toda lattice perturbed by a 
stochastic dynamics preserving energy and momentum. The strength of the stochastic noise 
is controlled by a parameter 7. We show that heat transport is anomalous, and that the 
thermal conductivity diverges with the length n of the chain according to «(n) ~ n a , with 
< a < 1/2. In particular, the ballistic heat conduction of the unperturbed Toda chain is 
destroyed. Besides, the exponent a of the divergence depends on 7. 

1 Introduction 

For a one dimensional system of length L, the thermal conductivity can be defined through the 
stationary flux of energy induced by connecting the system to two thermostats at different tem- 
peratures Ti and T r . This flux of energy Jl is proportional to the difference of temperature 
AT = Ti — T r , and we define the thermal conductivity kl as 

AT 

Jl = K L —. (1) 

If lim lim kt j = k exists and is finite, then the conductivity is normal and the system is said 

L^ooAT-0 

to satisfy Fourier's law [7]. The limit k is the thermal conductivity of the system. 

It is well known, by numerical experiments and certain analytical considerations, that the ther- 
mal conductivity diverges for one dimensional systems of oscillators with a momentum conserving 
dynamics [TS1 E] • This is also consistent with some experimental results on the length depen- 
dence of the thermal conductance of carbon nanotubes [231 E3 • In the case of a chain of harmonic 
oscillators, Jl can be computed explicitly [TI5] and does not decrease when the size L of the system 
increases. This is due to the ballistic transport of energy carried by the non-interacting phonons, 
and it happens also for optical (pinned) harmonic chains, where momentum is not conserved. Bal- 
listic transport of energy is expected for all systems whose dynamics is completely integrable [23] , 
as for the Toda chain [201 H5]. 

Numerical evidence shows that non-integr ability (whatever the definition of this concept one 
considers) is not a sufficient condition for normal conductivity, in particular for anharmonic chains 
of unpinned oscillators like the FPU model [14]. An energy superdiffusion is expected in these 
momentum conserving systems, and the thermal conductivity diverges as kl ~ L a , for some 
a 6 (0,1). There exists a wide debate in the physical literature about the existence or non- 
existence of one (or more) universal value of a [T71 [TS] . 
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Stochastic perturbations of the dynamics have been introduced in order to understand these 
phenomena. A model where Langevin thermostats are attached to each oscillator in a harmonic 
chain was first introduced by Bolsterli et al 0]. This noise destroys all types of conservation laws, 
including energy, and the corresponding conductivity is finite (see [5] and [5] for the anharmonic 
case). Also, stochastic perturbations that preserve only the energy of the system give finite thermal 
conductivity [3]. 

The situation changes dramatically if the stochastic perturbation conserves both the energy 
and the momentum pQ. Such stochastic conservative perturbations model the chaotic effect of 
non-linearities. These systems may then be seen as completely non-integrable, since the only 
conserved quantities left are the energy and the momentum. In this case, at least for harmonic 
interactions, the thermal conductivity can be explicitly computed by the Green-Kubo formula. 
For unpinned models, it remains finite only in dimension d > 3, while it diverges when d = 1 or 2. 
More precisely kl ~ \[L when d = 1, and kl ~ logL when d = 2 [H[2]. Analytical considerations 
for the same harmonic stochastic systems in non-equilibrium setting give the same results for the 
thermal conductivity computed according to (fT]), see 16J . For anharmonic interactions, rigorous 
upper bounds can be established, again by the Green-Kubo formula. In the one dimensional case, 
this leads to < C\J~L. 

These rigorous results motivated us to analyze the effect of stochastic perturbations on another 
completely integrable system, the Toda chain. The thermal conductivity of Toda lattices was 
already studied in [13] . where it was shown that the ballistic energy transport is destroyed for a 
diatomic system. In contrast with the harmonic case where many computations can be performed 
analytically, the nonlinear dynamics considered here has to be solved numerically. We considered 
a chain in the nonequilibrium steady state setting, with two Langevin thermostats at different 
temperatures attached to its boundaries. We chose the simplest possible stochastic perturbation 
conserving both momentum and energy: each couple of nearest neighbor particles exchange their 
momentum at random times distributed according to an exponential law of parameter 7 > 0. 

Our main results are the following: 

A) As soon as some noise is present, i.e. 7 > 0, the ballistic transport is immediately destroyed 
(as in the harmonic case) and energy superdiffuses, with kl ~ L a for < a < 1/2; 

B) the exponent a seems to depend on the noise strength 7, and is increasing with 7. 

If [Al was somehow expected, [B] is quite surprising. It may be explained by the noise destroying 
some diffusive phenomena due to non-linearities, like localized breathers, with the result that 
current-current correlation decays slower when more noise is present. Besides, [B] suggests that 
any theory claiming the existence of a universal parameter a has to be properly circumstanced. 

This paper is organized as follows. The dynamics we consider is described in Section [21 The 
numerical simulations we performed are described in details in Section [3[ whereas the obtained 
results are discussed in Section [U 

2 The stochastic dynamics 
2.1 Description of the system 

Hamiltonian. The configuration of the system is given by {q%,Pi, i — 0, . . . , n} 6 R 2n+2 , where 
qi is the displacement with respect to the equilibrium position of the i-th particle, and pi is its 
momentum. All masses are set equal to 1. The Hamiltonian is given by 




(2) 



where the interaction potential is defined by 



V(r) = T e br +ar + c. 
b 



(3) 
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and a,b,c are constants, a > 0, b > 0. Observe that if the product ab is kept constant, the 
harmonic chain is obtained in the limit b — > 0, while the hard sphere system is recovered as 
b — > oo. In our simulations, we chose a = 1/b and c = —1/b 2 in order for V to be non-negative 
and to be minimal at r = 0. The potential is therefore determined by a single parameter 6, which 
determines the strength of anharmonicity. 

Boundary conditions. We set qo = 0, which amounts to removing the center-of-mass motion 
by attaching the particle at the left end of the system to a wall. However, we do not fix the total 
length q n , and consider free boundary conditions on the right end. We checked that our numerical 
results are robust with respect to the boundary conditions. In particular, the same kind of scalings 
are obtained for fixed boundary conditions. 

2.2 Description of the dynamics 

The stochastic dynamics we consider has the following generator: 

L = A + + B n ) +-/S + rd Pn , (4) 
where £ and 7 are two positive constants. In 0$, A is the Hamiltonian vector field: 

n 

A = J2(Pid qi -d qi Hd Pi ) , (5) 
Bj are the generators of the Langevin thermostats attached at atom j = 1 and j = n: 

11 pa. ■ (6) 

and S is the generator of the random exchanges of momenta between nearest neighbor atoms: for 
any smooth function /, 

n-l 

Sf(q,p) = J2(f^p l ' %+1 )-.f(q,p)), (7) 

i=l 

where p l ' l+1 e R n is defined from p g E™ by 

Pi = Pi+i , Pi+i = Pi , (8) 

and p l j 1+1 — pj if j 7^ i, i + 1. Finally r is the strength of a constant external force applied to the 
last particle n. In ([S]), we choose the temperatures Tj = \ = T\ and Tj =n = T r . 

2.3 Energy currents 

We define the energy of the oscillator i for 1 < i < n — 1 as 

£ i = Y + \ " qi - l) + v{qi+1 ~ qt) ) ■ (9) 

Locally the energy conservation is expressed by the stochastic differential equation 

d£i (*) = dJi- hi {t) - dJi, i+1 (t) . (10) 

The energy currents J^i+i (t) are the sum of contributions from the Hamiltonian and the stochastic 
mechanisms. For i = 1, . . . , n — 2, the currents are 

Ji,i+i{t) = [ (j? ii+1 +7J? )<+ i) ds + Ml %+1 {t) , (11) 
Jo 
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where M^ i+1 (t) is a martingale, 

1 



JM+i = - J (P< +fi+i)V"(ft+i - ft) (12) 
is the instantaneous Hamiltonian current, while 

JM+i = | (rf-rf+i) (13) 

is the instantaneous stochastic current (the intrinsic transport of energy due to the stochas- 
tic exchange). The martingale term M~ t i+1 (t) can be characterized in the following way: let 
{iVi.i + i(i)}™r l 2 be independent Poisson processes of intensity 7. Then 

Ml l+1 {t) = J \(p 2 l { S -)-pi +1 { S -)){dN l , l+1 { S )~ 1 d S ) ) 

where p^(s~) = lim (i). At the boundaries of the system, the energy currents are 

Jb,i(*) = / |(r,-p?(s))d« + Vi^ / Pl(«)«M*), 



(I 



<-^n— l,n(^) 



I (Pn( s )- T r) +Tp„(s) 



ds + v£T r / Pn(s) dw n (s) , 



where u>i(£) and w n (t) are independent standard Wiener processes, and the last integrals on the 
right hand side of the previous formulas are Ito stochastic integrals. 

2.4 The stationary state 

If Ti = T r = T, we know explicitly the stationary probability measure of the dynamics, given by 
the Gibbs measure 

( — H+rq )/T n — B/T n—1 

^^rV II = i^T) n (e^'^MP.) dr ndPn , (14) 

where 7*j = qi — is the relative displacement, Z n (T,r) is a normalization constant, and 
£ = p, 2 j/2 + y(ri)/2 + V(r n )/2 — rr„ is a boundary term. 

If T; ^T r , there is no explicit expression of the stationary measure for anharmonic potentials. 
For certain classes of anharmonic potentials, the results of [TS1 El show that there exists a unique 
stationary probability measure. The assumptions on the potential made in [9] or similar works 
are not satisfied by the Toda potential (in particular, the growth at infinity is too slow in the limit 
r — > +00), but we believe that the techniques from [HI IS] can be extended to treat the case under 
consideration here. 

We denote by ( • ) the expectation with respect to this stationary measure, as well as the 
expectation on the path space of the dynamics in the stationary state. By stationarity we have 

(Ji, i+ i(t))=t (j£ i+1 +7J!,i+i ) ='-tJ n . 

Because of energy conservation, J„ does not depend on i, but only of the size n of the system. 
Consequently, 

J «- n _2^ {iu+l)+ 2 n-2 ■ (15j 
»=i 

In view of ([T]), the thermal conductivity can be defined by 

k„(T,t) = lim — — ^— . (16) 

y ' ; t^two Ti-T r v ' 
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It is clear from (fT5|) -([T5 1) that the direct contribution of the stochastic current to the conductivity 
is close to 7 and remains bounded in n. Hence only the first term of (|15p . namely the Hamiltonian 
current 



n-2 

rliam 



n-2 

roEOw). ( 17 ) 



n-2 

i=i 

can be responsible for a possible divergence of the conductivity. In the sequel, we hence consider 
the conductivity 

rham 

k% ub1 (T,t) = lim — 12 

T r — >T 1 ' 

rather than (fTtj|) . We are also motivated by the following numerical considerations. As reported 
in the sequel, we numerically observe that K^ am ~ n a for some a £ (0,1), hence J^ am ~ n v 
for v — a — 1 £ (—1,0). As a consequence, the second term of (fT5|) is indeed negligible with 
respect to the first term, in the limit n — * oo. The regime of large n may yet be difficult to reach 
numerically, so that the stochastic current contribution in (|15[) may be small but not negligible 
compared with (fT7|) for the considered values of n. 

By a linear response theory argument, the thermal conductivity of the finite system can also 
be defined by a Green-Kubo formula: 



-i />oo n ~2 

jrJ £GW*)j?,3(°)>Tr *• 



^ K (T,t) = ^I £<j?, i+1 (t)j?, 2 (0)> TT + (18) 



where here the expectation ( • ) T is with respect to the dynamics starting with the equilibrium 
Gibbs measure given by (fT^|) (we assume here that the integral fTS]) indeed exists). 

In principle, for finite n, 7^ k„, but we expect that they have, qualitatively, the same 
asymptotic behavior as n — > 00. 



3 Numerical simulations 

3.1 Implementation 

All the simulations performed in this work were done with XJ = 1.05 and T r = 0.95. This 
temperature difference is small enough so that the thermal conductivity around T = 1 should be 
approximated correctly. We also set the external force to r = 0. 



3.1.1 Integration of the dynamics 

We denote by g^iPf approximations of qi(t m ),pi(t m ) at time t m = mAt. The time-discretization 
of the dynamics with generator (U) is done with a standard splitting strategy, decomposing the 
generator as the sum of A + £(.Bi + B n ) and jS. 

The Hamiltonian part of the dynamics and the action of the thermostats on both ends of the 
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chain are taken care of by the so-called BBK discretization [5] of the Langevin dynamics: 



P 



m+l/2 



m+1 



m+1 



pT - ^ q Mi m ) 



+kn[-^Pn+\l^T r G™y 



At p. 



m+l/2 



P 



m+l/2 



At 



v 9s W i+1 ) 



r 



(19) 



where <5^i and 5i yTl are Kronecker symbols and G™, G™ are independent and identically distributed 
random Gaussian variables of mean and variance 1. Notice that the last step in the algorithm, 
written as an implicit update of the momenta, can in fact be rewritten in an explicit manner. Alter- 
natively, one can first integrate the Hamiltonian part of the dynamics with the Verlet scheme (21j 
using a time step At, and next analytically integrate the Ornstein-Uhlenbeck processes on the 
momenta at both ends of the chain, associated with the generator £(-Bi + B n ). In this work, we 
rather considered algorithm ([19)1 . We tested two different friction parameters, £ = 1 and £ = 0.1, 
to study how the results depend on the boundary conditions. 

The noise term with generator jS is simulated by exchanging pi and Pi+\ at exponentially 
distributed random times, with an average time 7 -1 between two such exchanges. More precisely, 
we attach to each spring a random time r™, with drawn from an exponential law with parameter 
7. This time is updated as follows: if t™ > At, then r" l+1 = t" 1 — At, otherwise pi and Pi+\ are 
exchanged and t™ +1 is resampled from an exponential law of parameter 7. 



3.1.2 Initial conditions and thermalization 

The initial conditions are chosen by imposing a linear temperature profile. We used to this end 
the Langevin dynamics of generator L\q = A + £B, where A is given by ([5]), and 

B = E T < - Pi d ^ . T i = ^-4*i + —T r . 

1=1 

Once a linear temperature profile, obtained from the time-average of the local kinetic energy, is 
indeed obtained, the term £B is switched off, and replaced by + B n ) + jS. The dynamics 
with generator (J3]) is then integrated using the numerical scheme described in Section [3.1. 11 and 
the spatially averaged instantaneous Hamiltonian current is monitored. At time t m = mAt, this 
current is defined as 

n-2 

r^Eta' (20) 



j 



i=l 



where the instantaneous Hamiltonian current J° j + i is defined as in l|12[). upon replacing qi(t) and 
Pi(t) by their approximations gf 1 and p™ . The thermalization time is somehow loosely defined 
as the time after which the variations of the instantaneous current stabilize (see Figure [1] for 
an illustration). This time could be determined more carefully by estimating some local-in-time 
variance of the current, and requiring that this variance stabilizes. 



G 




0.0 1.0e+06 2.0e+06 3.0e+06 

Time 

Figure 1: Instantaneous Hamiltonian current ([2TJ)) as a function of time, in the case 6 = 1,7=10 
£ = 0.1, At = 0.05, n = 2 16 = 65, 536. The thermalization time is T thm ~2x 10 6 . 



3.1.3 Computation of the energy currents 

Once some steady state has been reached, the instantaneous Hamiltonian current (120|) is computed 
at each time step, and an approximation of the Hamiltonian current (|17p is obtained as a time 
average: 

1 M 

^ = ^Ej m ( 21 ) 

m— 1 

In the limit of a large number of iterations M and for a small time-step At, we have J^f ~ </^ am - 
Using the current (|21|) . we define the conductivity /?„ and the exponent a by 

£M = ^n a (22) 

A priori, the exponent a depends on all the parameters of the model: the anharmonicity param- 
eter b, the magnitude £ of the coupling to the end thermostats, the noise strength 7, and the 
temperatures T) and T r . In our simulations, we have explored how the numerical results (and in 
particular the exponent a) depend on the first three parameters, and we have kept T; and T r fixed. 
We wish to point out that the results reported here already required an extremely large CPU time. 
Indeed, for the largest system considered (n = 2 17 ), several months were needed to integrate the 
dynamics with a small enough time step to ensure accuracy (At = 0.05 here), on a time T = 10 7 
long enough such that convergence of the time average (f2"Tj) is reached (see Section lB. 2.21 for more 
details). 

3.2 Error estimates 

There are two types of error in the numerical estimation of the currents: a systematic error (bias) 
due to the time step error (At > 0), and a statistical error due to the finiteness of the sampling 
(M < +00). We consider successively these two issues. 

3.2.1 Choice of time step 

The time step should be small enough in order for the dynamics to be numerically stable. When 
the size n of the system, the noise strength 7, the anharmonicity b or the Langevin friction £ are 
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increased, the time step should be reduced. Indeed, in all these cases, the energy of the system 
increases (at least locally). Due to nonlinearities, this energy may concentrate on a few sites, and 
hence trigger numerical instabilities. Such issues are not encountered with harmonic potentials, 
where some uniform stability condition is valid. 

In the case 6=1 and £ = 0.1, most of our computations have been done with At = 0.05. 
However, for the largest systems, and for the largest values of 7, we had to use the smaller value 
At = 0.025 (otherwise, the simulation blows up due to numerical instabilities, as for n = 2 14 and 
7 = 1). For 6=1 and £ = 1, that is a stronger noise at the boundaries thermostats, we also 
observe that we have to reduce the time step. We worked with At = 0.025 for all values of n. 
When 6 is increased from 6 = 1 to 6 = 10 (with £ = 0.1), the potential energy becomes stiffer, 
and we again need to use a smaller time step. In the case 6 = 10 and £ = 0.1, we worked with 
At = 0.01 for all values of n and 7, except for the large values of 7 and when n > 2 13 , for which 
we used At = 0.005. 

Let us now describe two artifacts of the numerical results that occur when the time step is 
too large. Observing them in practice is an indication that the time step is too large and should 
be reduced. In Figure [21 we plot the Hamiltonian current as a function of the chain length n, 
for 6 = 1, £ = 1 and 7=1, for two different time steps. These currents have been computed as 
the time averages (|21l) on simulations long enough. For At — 0.05, the current is not monotically 
decreasing with n, which is clearly a numerical artifact. The simulation blows up for n = 2 14 and 
is stable for n = 2 13 , but it is clear that the latter point cannot be trusted. 
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Figure 2: Hamiltonian current (|2"Tj) as a function of the chain length n, for At — 0.05 and 
At = 0.025 (with 6 = 1, £ = 1 and 7 = 1). 

On Figure [21 we plot the temperature profile, at the end of the simulation (T = MAt = 
4.24 x 10 7 ), in the case of a chain of length n = 2 13 (with 6 = 1, £ = 1 and 7 = 1). When we 
use At = 0.025, we obtain a decreasing temperature from the left end to the right end, which is 
in agreement with what is expected. When At — 0.05, the results are completely different, and 
physically unreasonable, which again shows that these results cannot be trusted. 

3.2.2 Variability of the results 

We mentioned earlier than we needed extremely long simulations to obtain a good accuracy. 
The reason for that can be well understood from Figure 01 on which we plot the instantaneous 
Hamiltonian current j m (defined by ((20)) ) and its time average (|2Tj) as a function of time, for a 
chain of length n — 2 (for the parameters 6 = 1, £ = 1 and 7 = 1). We observe that j m roughly 
oscillates between —0.02 and 0.02, whereas its time average is close to 10~ 4 . Hence the variability 
of j m is extremely large in comparison with its expectation. We hence need to run a very long 
simulation to be able to average out most of the fluctuation and obtain (j m ) with a good accuracy. 
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Figure 3: Averaged temperature profiles along the chain of length n = 2 13 = 8192, for At = 0.05 
and At = 0.025 (with b = 1, £ = 1 and 7 = 1) at the end of the simulation. 



These heuristic considerations can be quantified by a simple computation. The standard 
deviation of the instantaneous current in Figure S] is of the order of a ~ 0.02, while the average 
value of the current is /1 ~ 10~ 4 . In addition, it is possible to estimate the typical correlation time 
Tcorr using block averaging (also called batch means in the statistics literature), see [Til [12]. Here, 
we obtain r corr ~ 10 3 . The time i rcq required to obtain a 1% relative accuracy on the average 
current is such that 

= = 0.01 fi. 

y ^req/ "'con 

This yields t req ~ 4 x 10 11 . Since the time step is At ~ 0.05, this means that a huge number of 
time steps should be used to reach convergence. 

For large values of n and large values of 7, we observe that (j m ) is very small (see the numerical 
results below). Its accurate computation hence needs an even longer simulation time. In addition, 
the cost of the simulation of a chain, on a given time range, linearly increases with the size of the 
chain. This explains why computing the average currents for the longest chains is an extremely 
expensive task. 



instantaneous current — 
time-averaged current 
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Figure 4: Instantaneous Hamiltonian current (|2"0)) and its time average (|2ip as a function of time, 
for a chain of length n = 2 14 (for the parameters 6 = 1, £ = 1 and 7 = 1). 
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a 



a 



(6=1, £ = 1) (6 = 1, £ = 0.1) (6 = 10, £ = 0.1) 



0.001 

0.01 

0.1 

1 



0.10 
0.11 
0.32 
0.44 



0.17 
0.30 
0.44 



0.25 
0.32 
0.43 



Table 1: Conductivity exponent a, estimated from ([22)) . for different values of 7 and £, and 
different potential energies. 

3.3 Numerical results 

We present the conductivity ([2"2"|) as a function of the system length in Figure [5] We considered the 
choices 6=1 and b = 10 for £ = 0.1 to study the influence of the potential energy anharmonicity. 
We also simulated the system with 6=1 and £ = 1 to study how the results depend on £. The 
total simulation time for each point ranges from T = MAt = 10 6 to T — 5 x 10 7 , depending on 
the size n of the system. 

The slope a, defined by ([2"2")1 , is estimated using a least square fit in a log-log diagram. This 
estimate is quite sensitive to the choice of the number of points entering the fitting procedure, and 
only the very first digits of the estimated slope are reliable. Theoretical results (see [2]) show that 
the exponent a is expected to be lower than 0.5. Our numerical results, gathered in Table [TJ are 
in accordance with the theoretical upper bound. 

We also observe that, for 7 = 1, the value of a is close to 0.5. Now recall that, in the harmonic 
case V(r) = ar 2 , the value of a is always equal to 0.5, independently of 7. This seems consistent 
with the fact that, for large values of 7, the precise details of the potential V do not matter (the 
dynamics is mostly governed by the stochastic terms) , and the behaviour of the system is close to 
the harmonic behaviour. 

4 Discussion of the numerical results 

Several conclusions can be drawn from the numerical results given in the previous section: 

(i) The ballistic transport, observed in the deterministic Toda lattice, which is due to the com- 
plete integrability of the system [23] . is broken by the presence of noise in the dynamics. 
Energy transport becomes superdiffusive, i.e. n n ~ n a , with a € (0,0.5). For a low level 
of noise (7 small), this superdiffusive regime may be seen only for systems large enough 
(71 > 2 12 or more, depending on the stiffness parameter 6 of the system and the coupling £ 
to the boundaries thermostats). The asymptotic regime in n for the conductivity is attained 
for smaller values of n when 7 or £ is larger, or when 6 is smaller. 

(ii) The value of a seems to depend on the noise strength 7 in a monotonically increasing way. 
If had the same behavior, this would suggest that increasing the noise induces a slower 
time decay of the current-current correlations, in contradiction with the naive intuition that 
a stronger noise enhances the decay of time correlations. We believe that noise tends to 
suppress scattering effects due to the nonlinearity of the interaction. Observe that, in the 
harmonic case V(r) = ar 2 , the value of a is always equal to 0.5 and in particular does not 
depend on 7 [TJ [2J HE] . The dependence of the exponent a on the noise strength 7 in the 
Toda case also contradicts general theories about the universality of a (see for instance [12). 
These theories need then to be restricted to more specific dynamics. 

acknowledgements This work is supported in part by the MEGAS non-thematic program 
(Agence Nationale de la Recherche, France). 
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Figure 5: Conductivity K„ , defined by 1 ^, as a function of the syste 



m size n. 
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